El problema que vamos a tratar de resolver es predecir el salario de un jugador de la NBA para la temporada 2019-2020 en base a sus estadísticas de juego durante la temporada 2018-2019.
Las técnicas de regularización son útiles para trabajar con conjuntos con gran cantidad de variables, las cuales pueden introducir variabilidad en las estimaciones de los parámetros.
Los datos provienen de la página Basketball Reference y fueron previamente trabajados por nosotros para obtener el formato actual.
Las variables del set son:
diccionario = read_csv("diccionario_terminos.csv")
diccionario
En el glosario de Basketball Reference pueden encontrar una descripción más exhaustiva de cada una de estas métricas.
# Los datos de salario son para la temporada 2019-2020
nba = read_csv(" nba_player_stats_salary_2019_2020.csv") %>%
rename(salary = mean_salary_2019_2020) %>%
mutate(Pos = str_remove(string = Pos, pattern = "\\-.*")) %>%
mutate_all(~replace(., is.na(.), 0))
glimpse(nba)
## Rows: 395
## Columns: 50
## $ Player <chr> "Jaylen Adams", "Steven Adams", "Bam Adebayo", "LaMarcus Aldr…
## $ Pos <chr> "PG", "C", "C", "C", "SG", "C", "PF", "SF", "SF", "PF", "PF",…
## $ Age <dbl> 22, 25, 21, 33, 23, 20, 28, 25, 25, 30, 24, 34, 21, 24, 33, 3…
## $ Tm <chr> "ATL", "OKC", "MIA", "SAS", "UTA", "BRK", "POR", "ATL", "MEM"…
## $ G <dbl> 34, 80, 82, 81, 38, 80, 81, 48, 43, 25, 72, 10, 67, 81, 69, 8…
## $ GS <dbl> 1, 80, 28, 81, 2, 80, 81, 4, 40, 8, 72, 2, 6, 32, 69, 81, 70,…
## $ MP <dbl> 428, 2669, 1913, 2687, 416, 2096, 2292, 463, 1281, 322, 2358,…
## $ FG <dbl> 38, 481, 280, 684, 67, 335, 257, 64, 150, 21, 721, 49, 183, 1…
## $ FGA <dbl> 110, 809, 486, 1319, 178, 568, 593, 157, 276, 69, 1247, 121, …
## $ `FG%` <dbl> 0.345, 0.595, 0.576, 0.519, 0.376, 0.590, 0.433, 0.408, 0.543…
## $ `3P` <dbl> 25, 0, 3, 10, 32, 6, 96, 24, 9, 9, 52, 21, 67, 81, 145, 131, …
## $ `3PA` <dbl> 74, 2, 15, 42, 99, 45, 280, 77, 34, 40, 203, 64, 202, 217, 43…
## $ `3P%` <dbl> 0.338, 0.000, 0.200, 0.238, 0.323, 0.133, 0.343, 0.312, 0.265…
## $ `2P` <dbl> 13, 481, 277, 674, 35, 329, 161, 40, 141, 12, 669, 28, 116, 1…
## $ `2PA` <dbl> 36, 807, 471, 1277, 79, 523, 313, 80, 242, 29, 1044, 57, 202,…
## $ `2P%` <dbl> 0.361, 0.596, 0.588, 0.528, 0.443, 0.629, 0.514, 0.500, 0.583…
## $ `eFG%` <dbl> 0.459, 0.595, 0.579, 0.522, 0.466, 0.595, 0.514, 0.484, 0.560…
## $ FT <dbl> 7, 146, 166, 349, 45, 197, 150, 26, 37, 12, 500, 15, 36, 89, …
## $ FTA <dbl> 9, 292, 226, 412, 60, 278, 173, 35, 64, 16, 686, 22, 62, 102,…
## $ `FT%` <dbl> 0.778, 0.500, 0.735, 0.847, 0.750, 0.709, 0.867, 0.743, 0.578…
## $ ORB <dbl> 11, 391, 165, 251, 3, 191, 112, 24, 48, 18, 159, 9, 58, 27, 5…
## $ DRB <dbl> 49, 369, 432, 493, 20, 481, 498, 60, 203, 36, 739, 45, 139, 1…
## $ TRB <dbl> 60, 760, 597, 744, 23, 672, 610, 84, 251, 54, 898, 54, 197, 2…
## $ AST <dbl> 65, 124, 184, 194, 25, 110, 104, 23, 128, 19, 424, 5, 47, 269…
## $ STL <dbl> 14, 117, 71, 43, 6, 43, 68, 22, 54, 4, 92, 4, 46, 65, 91, 52,…
## $ BLK <dbl> 5, 76, 65, 107, 6, 120, 33, 13, 37, 1, 110, 7, 22, 4, 21, 4, …
## $ TOV <dbl> 28, 135, 121, 144, 33, 103, 72, 23, 58, 14, 268, 8, 55, 63, 1…
## $ PF <dbl> 45, 204, 203, 179, 47, 184, 143, 48, 112, 25, 232, 32, 140, 1…
## $ PTS <dbl> 108, 1108, 729, 1727, 211, 873, 760, 178, 346, 63, 1994, 134,…
## $ PER <dbl> 7.6, 18.5, 17.9, 22.9, 7.5, 18.5, 13.2, 11.2, 12.8, 4.6, 30.9…
## $ `TS%` <dbl> 0.474, 0.591, 0.623, 0.576, 0.516, 0.632, 0.568, 0.516, 0.569…
## $ `3PAr` <dbl> 0.673, 0.002, 0.031, 0.032, 0.556, 0.079, 0.472, 0.490, 0.123…
## $ FTr <dbl> 0.082, 0.361, 0.465, 0.312, 0.337, 0.489, 0.292, 0.223, 0.232…
## $ `ORB%` <dbl> 2.6, 14.7, 9.2, 10.3, 0.8, 9.6, 5.3, 5.3, 4.1, 6.1, 7.3, 3.3,…
## $ `DRB%` <dbl> 12.3, 14.8, 24.0, 19.8, 5.1, 24.0, 22.6, 13.9, 18.1, 12.5, 30…
## $ `TRB%` <dbl> 7.4, 14.7, 16.6, 15.1, 3.0, 16.8, 14.2, 9.5, 11.0, 9.2, 19.3,…
## $ `AST%` <dbl> 19.8, 6.6, 14.2, 11.6, 8.9, 7.9, 6.0, 6.9, 15.0, 7.8, 30.3, 2…
## $ `STL%` <dbl> 1.5, 2.0, 1.8, 0.8, 0.7, 1.0, 1.4, 2.2, 2.1, 0.6, 1.8, 0.7, 1…
## $ `BLK%` <dbl> 1.0, 2.4, 3.0, 3.4, 1.1, 4.5, 1.2, 2.4, 2.7, 0.3, 3.9, 2.0, 1…
## $ `TOV%` <dbl> 19.7, 12.6, 17.1, 8.8, 13.9, 13.0, 9.7, 11.8, 16.0, 15.5, 14.…
## $ `USG%` <dbl> 13.5, 16.4, 15.8, 26.9, 24.4, 15.9, 13.7, 17.2, 12.6, 12.0, 3…
## $ OWS <dbl> -0.1, 5.1, 3.4, 6.4, -0.4, 4.4, 3.0, 0.1, 0.8, -0.3, 8.9, 0.1…
## $ DWS <dbl> 0.2, 4.0, 3.4, 2.9, 0.4, 3.3, 2.8, 0.4, 1.8, 0.0, 5.5, 0.2, 1…
## $ WS <dbl> 0.1, 9.1, 6.8, 9.3, 0.0, 7.6, 5.8, 0.4, 2.7, -0.2, 14.4, 0.3,…
## $ `WS/48` <dbl> 0.011, 0.163, 0.171, 0.167, 0.002, 0.175, 0.121, 0.043, 0.100…
## $ OBPM <dbl> -3.8, 0.7, -0.4, 2.4, -4.2, 0.2, 0.1, -2.4, -1.6, -4.5, 6.3, …
## $ DBPM <dbl> -0.5, 0.4, 2.2, -0.6, -2.1, 1.4, 0.6, -0.3, 2.4, -1.9, 4.1, -…
## $ BPM <dbl> -4.3, 1.1, 1.8, 1.8, -6.3, 1.6, 0.7, -2.7, 0.8, -6.4, 10.4, -…
## $ VORP <dbl> -0.2, 2.1, 1.8, 2.6, -0.5, 1.9, 1.5, -0.1, 0.9, -0.4, 7.4, -0…
## $ salary <dbl> 131678, 25842697, 3454080, 26000000, 2429400, 2376840, 925800…
Existen 395 observaciones y 50 variables en el dataset.
Veamos como es la distribución de salarios según la posición en el juego. Se agregan las etiquetas de los jugadores que cobran mayores sueldos.
top_players = c("James Harden", "Stephen Curry", "Blake Griffin", "Chris Paul", "LeBron James", "Klay Thompson", "Jimmy Butler", "Gordon Hayward", "Kyle Lowry")
ggplot(nba, aes(Pos, salary, fill=Pos)) +
geom_boxplot() +
geom_text(aes(label=ifelse((salary>30500000) & Player %in% top_players,as.character(Player),'')),hjust=1.1,vjust=0, size=3) +
theme_bw() +
labs(title= "Boxplots: salarios y posición de juego", x="Posición", y="Salario")
Observamos que la distribución varía un poco entre posiciones y hay varios jugadores que son outliers según el criterio de Tukey.
Realizamos un correlograma entre todas las variables cuantitativas.
nba %>%
select_if(is.numeric) %>% # selección variables numéricas
ggcorr(., layout.exp = 5, hjust = 1, size = 3.5, nbreaks = 5, color = "grey50") + # graficamos correlacion pearson
labs(title='Correlograma de variables cuantitativas')
Observamos que existen relaciones de diversa magnitud y signo entre todas las variables.
Seleccionamos algunas variables y vemos sus relaciones usando ggpairs.
¿Cómo es la relación del salario con las demás variables?
nba %>%
select(salary, Age, PTS, GS, DRB, TRB, AST, BLK) %>%
ggpairs() + theme_bw()
Observamos que la correlación de todas estas variables con el salario es positiva pero de diversa magnitud.
También se ve que existen correlaciones muy fuertes entre ciertas variables. Por ejemplo, entre TRB (Total de Rebotes) y DRB (Rebotes Defensivos).
Vamos a probar un modelo lineal que incluya todas las variables (excepto al jugador y equipo) y obtener las estimaciones de los parámetros junto a su p-valor e intervalo de confianza.
Creamos el modelo y accedemos a información de los coeficientes e intervalos usando tidy.
# Eliminamos jugador y equipo
nba = nba %>% select(-c(Player, Tm))
# Modelo lineal
modelo_lineal = nba %>% lm(formula = salary~., data = .)
#Coeficientes
lineal_coef= modelo_lineal %>% tidy(conf.int=TRUE)
Graficamos los coeficientes estimados.
lineal_coef %>% filter(!is.na(estimate)) %>%
ggplot(., aes(term, estimate))+
geom_point(color = "forestgreen")+
geom_pointrange(aes(ymin = conf.low, ymax = conf.high), color = "forestgreen")+
geom_hline(yintercept = 0, lty = 4, color = "black") +
labs(title = "Coeficientes de la regresión lineal", x="", y="Estimación e Int. Confianza") +
theme_bw() +
theme(axis.text.x = element_text(angle=90))
Graficamos los p-valores de mayor a menor para evaluar la significatividad individual de los coeficientes estimados.
lineal_coef %>% filter(!is.na(estimate)) %>%
ggplot(., aes(reorder(term, -p.value), p.value, fill=p.value))+
geom_bar(stat = 'identity', aes(fill=p.value))+
geom_hline(yintercept = 0.05) +
labs(title = "P-valor de los regresores", x="", y="P-valor") +
theme_bw() +
theme(axis.text.x = element_text(angle=90)) +
scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )
Notamos que:
Hay ciertos coeficientes estimados que presentan una gran variabilidad pero la escala de las variables puede ocultarnos la verdadera variabilidad de los estimadores.
Pocas variables tienen coeficientes significativos: PosPG, TOV%, G y Age.
Existen cuatro variables cuyo coeficiente estimado es NA: 2P, 2PA, PTS y TRB.
lineal_coef %>% filter(is.na(estimate))
Cuando una variable se puede expresar como una combinación lineal de otra, el modelo lineal de R devuelve los valores de los coeficientes estimados como NA.
Para evitar los problemas que puede introducir la escala, reescalamos las variables con el comando scale y repetimos el ajuste para estos nuevos datos.
# Reescalamos las variables numericas
nba_scaled = nba %>% mutate_at(vars(-Pos), scale)
# Nuevo modelo lineal
modelo_lineal_scal = nba_scaled %>% lm(formula = salary~., data = .)
lineal_coef_scal = modelo_lineal_scal %>% tidy(conf.int=TRUE)
lineal_coef_scal %>% filter(!is.na(estimate)) %>%
ggplot(., aes(term, estimate))+
geom_point()+
geom_pointrange(aes(ymin = conf.low, ymax = conf.high), color = "forestgreen")+
geom_hline(yintercept = 0, lty = 4, color = "black") +
labs(title = "Coeficientes de la regresion lineal", subtitle="Variables escaladas", x="", y="Estimacion e Int. Confianza") +
theme_bw() +
theme(axis.text.x = element_text(angle=90))
lineal_coef_scal %>% filter(!is.na(estimate)) %>%
ggplot(., aes(reorder(term, -p.value), p.value, fill=p.value))+
geom_bar(stat = 'identity', aes(fill=p.value))+
geom_hline(yintercept = 0.05) +
labs(title = "P-valor de los regresores",subtitle="Variables escaladas", x="", y="P-valor") +
theme_bw() +
theme(axis.text.x = element_text(angle=90)) +
scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )
Ahora observamos que:
lineal_coef_scal %>% filter(is.na(estimate))
Obtemos la evaluación de ambos modelos ¿Cómo esperan que sean los valores de diagnóstico?
modelos_lineales = list(lineal = modelo_lineal, lineal_escalado = modelo_lineal_scal)
map_dfr(.x = modelos_lineales, .f = glance, .id="modelo") %>%
select(modelo, r.squared, adj.r.squared, p.value)
La alta cantidad de variables y la existencia de una alta correlación entre varias de ellas ocasionan que los coeficientes estimados tengan alta varianza y que muchos de ellos no sean significativos en términos estadísticos. Las técnicas de regularización pueden ayudarnos a mejorar esta situación.
La libreria glmnet nos permite trabajar con modelos ridge, lasso y elastic net. La función que vamos a utilizar es glmnet(). Es necesario que le pasemos un objeto matriz con los regresores y un vector con la variable a explicar (en este caso los salarios).
Fórmula:
\(ECM + \lambda{[\alpha\sum_{j}|\beta_j| +(1-\alpha)\sum_{j}\beta_j^2]}\)
\(\lambda\) controla toda la penalidad, mientras que \(\alpha\) controla la penalidad de la elastic net y tiende un puente entre Lasso y Ridge.
Con el parámetro \(\alpha\) indicamos con qué tipo de modelo deseamos trabajar:
Ridge: \(\alpha=0\)
Lasso: \(\alpha=1\)
Elastic Net: \(0<\alpha<1\)
Realizamos una partición entre dataset de entrenamiento y testeo usando la función initial_split del paquete rsample.
train_test <- nba %>% initial_split(prop = 0.7)
train <- training(train_test)
test <- testing(train_test)
En este caso vamos a trabajar con \(\alpha=1\).
¿Cuál es la penalización que introduce el modelo Lasso?
¿Cómo impacta esto en las variables?
# Vector con los salarios
nba_salary = train$salary
# Matriz con los regresores
nba_mtx = model.matrix(salary~., data = train)
# Modelo Lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=1, # Indicador del tipo de regularizacion
standardize = F) # Que esta haciendo este parametro?
lasso_coef = lasso.mod %>% tidy() %>% arrange(step)
lasso_coef
El comando plot nos permite realizar dos gráficos relevantes.
Gráfico de coeficientes en función del lambda
plot(lasso.mod, 'lambda')
Gráfico de coeficientes en función de la norma de penalización
plot(lasso.mod)
¿Qué muestra cada uno de estos gráficos?
Podemos realizar los gráficos para los valores de lambda también con ggplot.
g1=lasso_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line() +
theme_bw() +
theme(legend.position = 'none') +
labs(title="Lasso con Intercepto", y="Coeficientes")
g2=lasso_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line() +
theme_bw() +
theme(legend.position = 'none') +
labs(title="Lasso sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
Veamos un poco mejor aquellas variables que sobreviven para mayores valores de lambda ¿Qué tienen en común todas estas variables?
# Seleccionamos los terminos que sobreviven para valores altos de lambda
terminos_sobrevientes = lasso_coef %>%
filter(log(lambda)>16.5, term != "(Intercept)") %>%
select(term) %>% distinct() %>% pull()
# Graficamos
lasso_coef %>% filter(term %in% terminos_sobrevientes) %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line(size=1) +
geom_hline(yintercept = 0, linetype='dashed') +
theme_bw() +
labs(title="Lasso sin Intercepto", y="Coeficientes", subtitle= "\"Mejores\" variables") +
scale_color_brewer(palette = 'Set1')
Vemos que las variables que “sobreviven” para mayores valores de lambda son las que están medidas con una escala mayor.
glmnetExisten dos maneras de estandarizar las variables en glmnet.
Setear standardize = TRUE. Con esto se estandariza las regresoras y los coeficientes estimados estan en la escala original de la variable.
Pasar los conjuntos de datos estandarizados.
Replicamos todo el análisis previo pero para los datos estandarizados.
# Modelo lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=1, # Indicador del tipo de regularizacion
standardize = TRUE) # Estandarizamos
lasso_coef = lasso.mod %>% tidy() %>% arrange(step)
lasso_coef
Gráfico de coeficientes en función del lambda
plot(lasso.mod, 'lambda')
Gráfico de coeficientes en función de la norma de penalización
plot(lasso.mod)
Con ggplot
g1=lasso_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line() +
theme_bw() +
theme(legend.position = 'none') +
labs(title="Lasso con Intercepto", y="Coeficientes")
g2=lasso_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line() +
theme_bw() +
theme(legend.position = 'none') +
labs(title="Lasso sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
Veamos ahora cuáles variables sobreviven para mayores valores de lambda.
# Seleccionamos los terminos que sobreviven para valores altos de lambda
terminos_sobrevientes = lasso_coef %>% filter(log(lambda)>13.1, term != "(Intercept)") %>%
select(term) %>% distinct() %>% pull()
# Graficamos
lasso_coef %>% filter(term %in% terminos_sobrevientes) %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
geom_line(size=1) +
geom_hline(yintercept = 0, linetype='dashed') +
theme_bw() +
labs(title="Lasso sin Intercepto", y="Coeficientes", subtitle= "\"Mejores\" variables") +
scale_color_brewer(palette = 'Set1')
Observamos que ahora tenemos otro set de “mejores” variables.
¿Podemos decidir cuál es el valor óptimo de lambda?
Para elegir el valor óptimo de lambda, lo común es realizar cross-validation. La función cv.glmnet nos permite realizar esto de manera sencilla.
Al igual que para la función glmnet cuenta con los parámetros:
x: matriz de variables
y: vector de la variable a predecir
alpha: tipo de modelo
standardize: flag lógico para estandarizar las variables
Nuevo parámetro:
Salida Base
lasso_cv = cv.glmnet(x=nba_mtx,y=nba_salary,alpha=1, standardize = T)
summary(lasso_cv)
## Length Class Mode
## lambda 100 -none- numeric
## cvm 100 -none- numeric
## cvsd 100 -none- numeric
## cvup 100 -none- numeric
## cvlo 100 -none- numeric
## nzero 100 -none- numeric
## call 5 -none- call
## name 1 -none- character
## glmnet.fit 12 elnet list
## lambda.min 1 -none- numeric
## lambda.1se 1 -none- numeric
Brinda mucha información:
lambda: valor de lambda
cvm (Cross-validation mean): es la media del MSE (error)
cvsd (Cross-validation Standard Error): desvio estandar del MSE (error)
cvup y cvlo: Limite superior e inferior
nzero: Coeficientes distintos de cero
lambda.min: lambda para el cual el MSE (error) es mínimo
lambda.1se: lambda que se encuentra a 1 desvío estandar del lambda.min
glm.fit: incluye cantidad de variables, el valor de lambda y el porcentaje de deviance explicada por el modelo.
Si imprimimos el objeto tenemos:
lasso_cv
##
## Call: cv.glmnet(x = nba_mtx, y = nba_salary, alpha = 1, standardize = T)
##
## Measure: Mean-Squared Error
##
## Lambda Measure SE Nonzero
## min 168429 3.626e+13 3.368e+12 17
## 1se 1188231 3.930e+13 3.814e+12 6
Gráfico Base
plot(lasso_cv)
El gráfico nos muestra la media del MSE con su límite superior e inferior y la cantidad de variables que sobreviven para cada valor de lambda.
Usando Broom
Obtenemos la información del objeto lasso_cv con las funciones tidy y glance.
# Información de CV en dataframe con tidy
lasso_cv %>% tidy()
# Lambda minimo y lambda a 1 desvio estandar
lasso_cv %>% glance()
Seleccionamos el lambda óptimo para crear el modelo final.
# Selección lambda óptimo
lasso_lambda_opt = lasso_cv$lambda.min
# Entrenamiento modelo óptimo
lasso_opt = glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=1, # Indicador del tipo de regularizacion
standardize = TRUE, # Estandarizamos
lambda = lasso_lambda_opt)
# Salida estandar
lasso_opt
##
## Call: glmnet(x = nba_mtx, y = nba_salary, alpha = 1, lambda = lasso_lambda_opt, standardize = TRUE)
##
## Df %Dev Lambda
## 1 17 65.9 168400
# Tidy
lasso_opt %>% tidy()
# Glance (no es muy informativo)
lasso_opt %>% glance()
Han quedado 17 variables y el modelo explica el 65,9% de la deviance.
En este caso vamos a trabajar con \(\alpha=0\). Vamos a replicar lo que ya realizamos para Lasso.
¿Cuál es la penalización que introduce el modelo Ridge?
¿Cómo impacta esto en las variables?
En este caso, ya seteamos parámetro para estandarizar las variables. Si no lo fijasemos, el default sería de igual modo standarize=TRUE.
#Modelo ridge
ridge.mod=glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=0, # Indicador del tipo de regularizacion
standardize = TRUE)
#Coeficientes tidy
ridge_coef= ridge.mod %>% tidy()
ridge_coef
¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso?
Gráfico de coeficientes en función del lambda
plot(ridge.mod, 'lambda')
¿Qué ven de distinto en este gráfico respecto al que obtuvimos con la regresión Lasso?
Gráfico de coeficientes en función de la norma de penalización
plot(ridge.mod)
g1=ridge_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Ridge con Intercepto", y="Coeficientes")
g2=ridge_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Ridge sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
#cross-validation
ridge_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0, standardize = T)
Gráfico Base
plot(ridge_cv)
Seleccionamos el lambda óptimo para crear el modelo final.
# Selección lambda óptimo
ridge_lambda_opt = ridge_cv$lambda.min
# Entrenamiento modelo óptimo
ridge_opt = glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=0, # Indicador del tipo de regularizacion
standardize = TRUE, # Estandarizamos
lambda = ridge_lambda_opt)
# Salida estandar
ridge_opt
##
## Call: glmnet(x = nba_mtx, y = nba_salary, alpha = 0, lambda = ridge_lambda_opt, standardize = TRUE)
##
## Df %Dev Lambda
## 1 50 65.91 2810000
# Tidy
ridge_opt %>% tidy() %>% mutate(estimate = round(estimate, 4))
En este caso conserva todas las variables, obteniendo un porcentaje de deviance explicado muy similar a Lasso: 65,9%.
El modelo Elastic Net incorpora los dos tipos de penalización: Lasso (Norma L1) y Ridge (Norma L2). Define un compromiso entre las penalidades de lasso y ridge. El parámetro \(\alpha\) regula la importancia de cada penalización, cuanto más cerca de cero será más importante la penalización del tipo Ridge y más cerca de 1, la tipo Lasso.
En este caso vamos a trabajar con \(\alpha=0.5\). Vamos a replicar lo que ya realizamos para Lasso y Ridge, individualmente.
# Modelo elastic net
elastic.mod=glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=0.5, # Indicador del tipo de regularizacion
standardize = TRUE)
# Coeficientes del modelo
elastic_coef= elastic.mod %>% tidy() %>% mutate(estimate = round(estimate, 4)) %>% arrange(step)
elastic_coef
¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso y Ridge?
Gráfico de coeficientes en función del lambda
plot(elastic.mod, 'lambda')
¿Qué ven en este gráfico de distinto a los dos anteriores?
g1=elastic_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Elastic Net con Intercepto", y="Coeficientes")
g2=elastic_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Elastic Net sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
elastic_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0.5, standardize = T)
Grafico Base
Presten especial atención al eje superior ¿Qué está sucediendo?
plot(elastic_cv)
Seleccionamos el lambda óptimo para crear el modelo final.
# Selección lambda óptimo
elastic_lambda_opt = elastic_cv$lambda.min
# Entrenamiento modelo óptimo
elastic_opt = glmnet(x=nba_mtx, # Matriz de regresores
y=nba_salary, #Vector de la variable a predecir
alpha=0.5, # Indicador del tipo de regularizacion
standardize = TRUE, # Estandarizamos
lambda = elastic_lambda_opt)
# Salida estandar
elastic_opt
##
## Call: glmnet(x = nba_mtx, y = nba_salary, alpha = 0.5, lambda = elastic_lambda_opt, standardize = TRUE)
##
## Df %Dev Lambda
## 1 19 65.37 405700
# Tidy
elastic_opt %>% tidy() %>% mutate(estimate = round(estimate, 4))
Vamos a comparar la relación entre el porcentaje de deviance explicada y lambda para los tres tipos de modelos que realizamos.
ridge_dev = ridge_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
ggplot(., aes(log(lambda), dev.ratio)) +
geom_point() +
geom_line() +
geom_vline(xintercept = log(ridge_lambda_opt), color='steelblue', size=1.5) +
labs(title='Ridge: Deviance') +
theme_bw()
lasso_dev = lasso_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
ggplot(., aes(log(lambda), dev.ratio)) +
geom_point() +
geom_line() +
geom_vline(xintercept = log(lasso_lambda_opt), color='firebrick', size=1.5) +
labs(title='Lasso: Deviance') +
theme_bw()
elastic_dev = elastic_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
ggplot(., aes(log(lambda), dev.ratio)) +
geom_point() +
geom_line() +
geom_vline(xintercept = log(elastic_lambda_opt), color='forestgreen', size=1.5) +
labs(title='Elastic Net: Deviance') +
theme_bw()
plot_grid(ridge_dev, lasso_dev, elastic_dev)
Con los modelos óptimos que encontramos se puede calcular cuál es el RMSE en los datasets de training y testing, para decidir cuál es el modelo que minimiza el error en las predicciones.
# Definimos una función augment que funcione con glmnet
augment_glmnet = function(df,y, model) {
# Matriz con los regresores
formula = as.formula(str_c(y, "~."))
data_matrix = model.matrix(formula, data = df)
predictions = predict(model, data_matrix)
pred_colname = str_c("predicted",y, sep="_")
df[pred_colname] = predictions
return(df)
}
Vamos a agregar dos modelos para comparar: el salario promedio del set de entrenamiento y el modelo lineal múltiple clásico.
# Salario promedio del set de entrenamiento
salario_promedio = mean(train$salary)
# Modelo lineal
modelo_lineal = lm(salary~., data = train)
Realizamos las predicciones ambos modelos.
# Salario promedio
prediccion_modelo_nulo = tibble(salary = train$salary, predicted_salary = salario_promedio)
# Predicciones del modelo lineal
prediccion_modelo_lineal = augment(modelo_lineal) %>% mutate(predicted_salary = .fitted) %>% select(salary, predicted_salary)
Realizamos las predicciones de los modelos de glmnet.
# Lista de modelos
modelos_glmnet = list(lasso=lasso_opt, ridge=ridge_opt, elastic=elastic_opt)
# Predicciones de los modelos glmnet
lista_predicciones_training = map(.x = modelos_glmnet, .f = augment_glmnet, df=train, y="salary")
# Agregamos las otras predicciones a la lista
lista_predicciones_training = lista_predicciones_training %>% prepend(list(nulo=prediccion_modelo_nulo, lineal = prediccion_modelo_lineal))
Obtenemos el RMSE para el set de entrenamiento.
map_dfr(.x = lista_predicciones_training, .f = rmse, truth="salary", estimate="predicted_salary", .id="modelo")
¿Cuál es el modelo que realiza la mejor predicción? ¿Qué esperan que suceda con el RMSE en el set de evaluación?
Obtenemos el RMSE para los modelos en el set de evaluación
# Predicción del promedio
prediccion_modelo_nulo = tibble(salary = test$salary, predicted_salary = salario_promedio)
# Predicción modelo lineal
prediccion_modelo_lineal = augment(modelo_lineal, newdata = test) %>% mutate(predicted_salary = .fitted)
# Predicciones glmnet
lista_predicciones_test = map(.x = modelos_glmnet, .f = augment_glmnet, df=test, y="salary")
# Lista completa de predicciones
lista_predicciones_test = lista_predicciones_test %>% prepend(list(nulo=prediccion_modelo_nulo, lineal = prediccion_modelo_lineal))
# RMSE en el set de evaluación
map_dfr(lista_predicciones_test, rmse, truth="salary", estimate="predicted_salary", .id="modelo")
¿Cuál es el modelo de mejor performance en el set de evaluación? ¿Qué sucedió con el modelo de regresión clásico?